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Order parameters based on spherical harmonics and Fourier coefficients aheady play a signifi- 
cant role in condensed matter research in the context of systems of spherical or point particles. 
Here, we extend these types of order parameter to more complex shapes, such as those encountered 
in nanoscale self-assembly applications. To do so, we build on a powerful set of techniques that 
originate in the computer science field of "shape matching." We demonstrate how shape match- 
ing techniques can be applied to identify unknown structures and create highly-specialized ad hoc 
order parameters. Additionally, we investigate the special symmetry properties of harmonic de- 
scriptors, and demonstrate how they can be exploited to provide optimal solutions to certain classes 
of problems. Our techniques can be applied to particle systems in general, both simulated and 
experimental, provided the particle positions are known. 



Quantitative measures of symmetry and order, such as 
order parameters and correlation functions, are often ap- 
plied within the physical and chemical sciences to study 
the structural properties of particle systems. Structural 
quantities are particularly important in condensed sys- 
tems, including nano and colloidal scale self-assembly ap- 
plications, where subtle differences in particle ordering 
can greatly affect the thermodynamical, physical, chem- 
ical, electrical, and optical properties of a material or 
device pitZi- Some particularly useful order parameters, 
known as "bond order parameters," were introduced by 
Nelson and co-workers in the context of 2d and 3d sim- 
ulations of point particles [HI, fQ]. These order parameters 
have since been widely applied to both simulated and ex- 
perimental systems for quantifying crystal-like ordering 
in systems of spherical particles. Some common appli- 
cations of bond order parameters include, but are not 
limited to, identifying small ordered clusters ^9till^, con- 
structing static[12-17 and temporal [181 HI] correlation 
functions, identifying structural defects [20|, and studying 
nucleation and growth [151 [16]. Bond order parameters 
have the advantage that they can give a representative 
value for a given structure regardless of spatial orienta- 
tion. Additionally, they are robust under random per- 
turbations due to thermal noise and highlight important 
rotational symmetries. 

Since bond order parameters were originally designed 
to quantify order in small point clusters [8, 9 , they can- 
not be applied directly to complex structures or parti- 
cle shapes. Thus, they fail, in many cases, to fully de- 
scribe the complex structures that arise in contempo- 
rary disciplines, such as nanoscale and colloidal assem- 
bly, soft matter physics, and the biological sciences. The 
field of nanoparticle assembly in particular encompasses 
a vast range of the structural complexity that is possi- 
ble for particle systems [6l [71 EH [22]. Here, nanometer 
and micron-sized colloidal particles with a wide range 
of shapes, compositions and interparticle forces self- 
assemble into unique structures, such as complex crys- 
tals reminiscent of atomistic condensed matter [23H26] . 
phase-separating structures similar to those observed for 



block copolymer and surfactant svstems [2T| [26H3T] , and 
hierarchical assemblies that resemble certain biological 
structures j32l [33] . Because there are no standard ac- 
cepted structural metrics or order parameters for describ- 
ing these systems, ad hoc analyses or visual inspection are 
often employed instead. This often yields relatively inac- 
curate or incomplete results when compared to statistical 
analysis. 

Here, we generalize and extend bond order parameters 
for applications involving complex structures. Our pri- 
mary focus is on the field of self-assembly, where we pre- 
dict that these new structural characterization schemes 
are immediately required and readily applicable. How- 
ever, the structural metrics that we introduce, known 
as "harmonic descriptors," can be applied across diverse 
fields, including soft matter, macromolecular and bio- 
logical sciences, and other applications where complex 
particle structures are encountered. To derive general 
structural metrics for these systems, we draw strongly 
from the computer science field of shape matching [34l [35] , 
which aids us both in constructing mathematical rep- 
resentations for describing complex structures, and in 
quantifying structural similarity based on these repre- 
sentations. Our approach involves constructing "shape 
descriptors" based on a collection of many bond order 
parameters computed over a range of length scales. By 
quantifying the similarity between pairs of these shape 
descriptors, we can derive order parameters and correla- 
tion functions that are applicable to complex structures. 
In this respect, the order parameters that we introduce 
actually represent the degree to which two collections of 
"sub-order parameters" match. This pairwise compara- 
tive approach to constructing order parameters is neces- 
sary for applications involving complex structures, which, 
unlike simple clusters of spherical particles, require a 
range of structural metrics to describe them completely. 
Shape matching techniques based on the harmonic de- 
scriptors described here have already been applied to 
complex particle systems in the context of fast database 
searches for retrieving macromolecules and proteins [36]- 
^38^. We have applied similar database searches in the 
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context of characterizing local structure in nanoparticle 
assemblies [1 13 [23]. In addition to this type of applica- 
tion, we demonstrate how harmonic descriptors can be 
applied to the broad range of structural characterization 
problems to which bond order parameters have tradition- 
ally been applied for simpler systems. 

This article is organized as follows. In section [Tj we 
describe how to extract patterns from complex particle 
systems that can be described by harmonic shape descrip- 
tors. In section |ll| we describe, in detail, how the descrip- 
tors can be computed mathematically, and explore their 
unique properties, such as rotational invariance and sen- 
sitivity to rotational symmetries. In section |IV[ we de- 
scribe how the harmonic descriptors can be used to solve 
representative problems from the field of computational 
and experimental self-assembly and computational biol- 
ogy. In section |V| we explore how the special symmetry 
properties of harmonic descriptors can be applied to solve 
unique problems that are not easily solved by other types 
of shape descriptors. The methods that we describe here 
are applicable to all types of particle systems, both sim- 
ulated and experimental, for which particle positions are 
known or can be accurately imaged. To aid with the dis- 
semination of these techniques and new algorithms using 
harmonic shape descriptors as a basis, we provide a soft- 
ware library via the web [39 . A general review of the 
topic of applying shape matching algorithms to particle 
systems can be found in References [34] and [35] . 



I. PATTERN EXTRACTION 

The first step in the structure characterization schemes 
that we will employ is to extract a representative struc- 
tural pattern from the particle system that can be "in- 
dexed" (i.e., described mathematically) by a harmonic 
shape descriptor. While this is relatively trivial for small 
clusters of spherical point particles [9 , for complex struc- 
tures, some physical and mathematical intuition is often 
required. As we will demonstrate in detail in section [ll| 
the harmonic shape descriptors that we introduce are 
best suited to index patterns on the unit circle, sphere, 
disk or ball. Such representations are often sufficient to 
distinguish between even very similar structures with a 
high degree of precision. The patterns may represent the 
particles themselves or some interesting pattern formed 
by their positions or density profile. In some cases, the 
raw data may be preprocessed to better extract impor- 
tant structural features, for example, by spatial coarse- 
graining, time averaging, or potential energy minimiza- 
tion. 

Structural patterns in general can be described by a 
set of positions {x} = {xi, X2, . . . x^} and corresponding 
weights {/} = {/i, /2, . . . /n}- For point cloud data, or 
raw particle coordinates, {x} represents the particle po- 
sitions and the weights {/} are equivalent and usually 
taken to be 1. For voxel data (i.e., volumetric data, of- 
ten used to describe density maps), {x} represents the 
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FIG. 1: Extracting local patterns, (a) Clusters of point par- 
ticles are expressed as patterns on the perimeter of the circle 
or on the surface of the sphere. (6) More complex particles 
are often idealized as point particles to the same effect. (6) 
Point clusters with r-dependence or spatial density maps can 
be decomposed into independent patterns for multiple radial 
shells. 



positions of bins on a grid with weights {/}. The same 
representation can be used to describe experimental im- 
ages; in this case, {x} represents the positions of the 
pixels and {/} represents their intensity. This notation 
allows us to write general equations for shape descriptors 
in section ITIl 



A. Local Structures 

As originally shown in the context of bond order 
parameters'^, a cluster of point particles is one simple 
type of structure that can be trivially represented by the 
projection of the points onto the surface of a circle or 
sphere. The patterns for two different point clusters are 
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shown in Fig. [T^. Notice that we exclude the center par- 
ticle from the pattern, since it has no specific direction 
relative to itself. The clusters need not strictly consist of 
point particles; so long as particle shape is not important, 
the structural pattern can be described by placing points 
at the particle centroids, as depicted in Fig. [T]3. Often, 
local structures are isolated from the bulk system by ap- 
plying a clustering algorithm. One standard scheme is 
to cluster all particles within a cutoff range [Q) [TT]. More 
specialized schemes can be applied for specific applica- 
tions, such as particles with complex shape [23 . 

Projecting patterns onto the circle or sphere neglects 
radial information. Therefore, such projections can lead 
to non-distinguishing patterns for structures with radial 
dependence. One solution is to decompose structures into 
concentric shells, projecting each shell onto the circle or 
sphere independently, as depicted in Fig. [l]^. Structures 
can be compared by matching corresponding shells. Al- 
ternatively, complex patterns can be represented on the 
unit disk or unit ball. These representations account for 
radial dependence as well as angular dependence. As de- 



scribed in section II B , representation on the disk or ball 
is optimal when decomposition into shells is inaccurate, 
or gives degenerate representations. 



B. Global Structures 

Structures with long-range ordering, such as crystals 
and phase- separated structures, cannot be directly in- 
dexed on the circle, sphere, disk or ball. Rather global 
patterns are constructed by combining different pieces 
of local information. One such method, originally in- 
troduced in the context of bond order parameters [9 , is 
based on computing the superposition of all local pat- 
terns in a sample [42 . Fig. [2|l shows the superposition 
of local patterns in a face centered cubic crystal, where 
local patterns are defined by the projection of neighbor 
directions on the unit sphere. Since crystals have long- 
range orientational ordering, the neighbor directions are 
coincident throughout the sample. This type of struc- 
tural pattern is independent of the shape of the underly- 
ing particles, and thus can be applied to many assembled 
structures, such as crystals of patchy particles [43], poly- 
hedral particles [23], or phase-separated structures that 
form micelles or cylinders arranged in crystalline super- 
lattices, such as dendrimers[44j. block copolymers [3 , or 
tethered nanoparticles [TT| |4Q| l45ti47] . This superposi- 
tion scheme is also applicable to many non-crystalline 
global phase-separated structures, such as layered and 
network structures formed by block copolymers or teth- 
ered nanoparticle svstems[4Ql. \4E\ l46] . Fig. [2]3 shows the 
superposition of local density maps represented on the 
unit ball for a phase-separated sheet structure formed by 
tethered nano spheres [46 . 

For structures without long-range orientational order- 
ing, the superposition of local patterns results in a ran- 
dom pattern over long ranges. Fig.[2]3 shows the superpo- 



sition of local patterns for an atomic liquid, which results 
in a uniform distribution on the sphere. Since the same 
pattern is inherent to all liquids, gases, gels, etc., regard- 
less of the underlying particle shape, superposition gives 
no structural information, other than indicating an ab- 
sence of long-range orientational order. Orientationally- 
disordered structures of all types can be differentiated 
by considering the probability distribution functions of 
local patterns, rather than the superposition. The prob- 
ability distribution scheme is also useful when compar- 
ing structures with complex local patterns where super- 
position becomes degenerate or non-distinguishing, such 
as complex network structures [411. This is depicted for 
the double-gyroid structure formed by tethered nano- 
rodsl401|41i in Fig.|2]i. 



II. HARMONIC DESCRIPTORS 

Given a pattern on the unit circle, sphere, disk or ball, 
the harmonic shape descriptors reviewed in this section 
can be used to index the pattern into a compact vec- 
tor representation. This vector, or "shape descriptor," 
can be compared with other shape descriptors to ob- 
tain a quantitative measure of similarity between struc- 
tures. In the derivation that follows, we introduce har- 
monic descriptors from a perspective that draws on el- 
ements of shape matching and signal processing. This 
contrasts with the physics-based perspective presented 
in the original derivation of bond order parameters [9j. 
The alternate perspective is meant to highlight the fun- 
damentally different way in which the structural metrics 
are used; whereas bond order parameters are often ap- 
plied directly as order parameters, harmonic descriptors 
represent structural fingerprints that must be matched to 
obtain structural information, which can subsequently be 
used for constructing order parameters. 

Before we introduce harmonic descriptors, it is impor- 
tant to understand why they are more useful than simpler 
shape descriptor methods. Consider, for example, the 
problem of mathematically comparing two simple struc- 
tures, such as the clusters shown in Fig. [l^. Perhaps 
the most obvious way to describe the different struc- 
tural patterns is to simply use the coordinates them- 
selves as the shape descriptor. While this greatly sim- 
plifies the initial step of creating a shape descriptor, it 
complicates matching significantly, since we do not typi- 
cally know a priori the optimal correspondence between 
the coordinates in different lists. Reordering the lists is 
an optimization problem that can be solved, for exam- 
ple, by applying the Hungarian method [48 . This type 
of problem scales as 0{N^) and thus quickly becomes 
inefficient for large N. A simple solution to the as- 
signment problem is to create a probability histogram 
on the unit circle, sphere, disk or ball, as depicted in 
Figs. [1] and [2j Since the histogram bins are indepen- 
dent of the order of the particles in the list (and number 
of particles), no assignment is required. Although this 
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FIG. 2: Extracting global patterns, (a) Superposition of local patterns for an fee erystal. (6) Superposition of local patterns 
for a phase separated lamellar structure formed by tethered nanospheres 11 . (c) Superposition and probability distribution of 
local patterns for a liquid. Since the structure has no long range orientational ordering, superposition is non-distinguishing, (d) 
Superposition and probability distribution of local patterns for a double-gyroid structure formed by tethered nano rods j40ll4T] . 
For this complex structure, the superposition of local patterns is non distinguishing, even though the structure has long-range 
ordering. 



representation, known as the "shape histogram [49^," is 
very useful for some applications, one drawback is that 
to compare patterns in a way that is rot at ion- invariant, 
the patterns (or the histograms) must be aligned prior 
to matching. This "registration f50j EJ" step is com- 
putationally expensive and potentially inaccurate if ap- 
plied naively. One elegant solution to these problems 
is to compute the discrete Fourier transform (DFT) for 
each shell in the shape histogram [52 . The DFT trans- 
forms the pattern into its frequency-domain representa- 
tion, which can be used to obtain rotation-independent 
harmonic descriptors through a simple mathematical op- 
eration as described in the following section. As an ad- 
ditional advantage, harmonic descriptors have adjustable 
frequency parameters that can be tuned to highlight im- 
portant rotational symmetries or give a variable degree 
of coarse-graining. The manner in which we compute the 
harmonic descriptors depends on the coordinate system 
best used to describe the structure. In the following sec- 
tions, we first introduce Fourier descriptors |52H54] . which 
are suited for indexing shapes on the unit circle {0 depen- 
dence) or sphere (O^cj) dependence), and then introduce 



Zernike moments which are suited for indexing shapes 
on the unit disk (r, dependence) or ball (r, 6>, (j) depen- 
dence). 

A. Fourier Descriptors 

Fourier descriptors are designed to efficiently index 
structural patterns on the circle or sphere. A Fourier 
transform represents a function as a sum of harmonic 
components. For a pattern along the Id perimeter of 
the circle, this representation is performed by summing 
complex exponential terms: 

fi^j) =^^i exp [iWj] j = 1, 2, ...n,,n. (1) 

i=0 

Here, f{Oj) is the intensity of the pattern at a particular 
point along the perimeter of the circle Oj. The terms tpi, 
known as "Fourier coefficients," indicate the strength of 
the pattern for a particular frequency i. Typically, we 
cut off the frequency i at some finite value imax^ since 
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FIG. 3: Fourier descriptors, (a) Decomposition of the angular 
pattern obtained for a 2d hexagonal cluster into a sum of sines 
and cosines. Since the cluster has 6-fold symmetry, the lead- 
ing coefficients are given by multiples of ^ = 6. (6) Decompo- 
sition of the angular pattern obtained for a 3d fee cluster into 
a sum of spherical harmonics. Since the cluster has 4 and 6- 
fold symmetry, the leading coefficients are given by multiples 
of ^ = 4 and i = 6. (c) Schematic of the resulting Fourier 
descriptors. The value of the each component in the descrip- 
tor is proportional to the contribution of the corresponding 
harmonic to the overall sum. Notice that rotation-invariant 
Fourier descriptors contain real, positive components, while 
rotation-dependent descriptors contain complex components. 



the information for high-frequency i becomes increas- 
ingly dominated by noise in the structure. The i = 
and £ = 1 terms only contain information regarding the 
position and center of mass of the pattern and are some- 
times excluded. If the pattern consists of more than one 
radial shell, we compute the Fourier transform for each 
shell independently. 

The Fourier coefficients contain structural information 
that can be used to create shape descriptors with various 



properties. The coefficients are given by: 



= — ^ fj exp [uejY 



(2) 



The term 6j is the angle of an input point Xj with in- 
tensity fj from our pattern {x},{/}. The coefficients 
ijji are complex numbers. The * denotes the complex 
conjugate. As outlined in the previous section, the input 
points from our pattern can represent either the positions 
of bins in the circular shape histogram or the raw data 
if no binning is performed. The representations become 
equivalent as the bin size approaches zero and only one 
point can occupy a given bin. 

Since the Fourier transform is a frequency-domain rep- 
resentation of the pattern, it gives the strongest signal 
for frequencies that reflect periodicities in the pattern 
around the circle. That is, patterns with n-fold rota- 
tional symmetry yield high values of ifji (and, as we will 
discuss later, q^). Since rotational symmetries are rela- 
tively insensitive to small changes to the pattern, Fourier 
descriptors are relatively insensitive to thermal noise, 
particularly for low- frequency coefficients. Although the 
Fourier coefficients in their complex number form are not 
rotation-invariant, we can convert them to their invari- 
ant form by computing the magnitude of each coefficient. 
The invariant circular coefficients are given by: 



1^,1 = ^^ri = ^MW+Wef 



(3) 



The Fourier invariants are positive real numbers. 

To obtain some physical understanding of the proper- 
ties of Fourier coefficients, consider the small cluster in 
Fig. [4^. If we normalize the cluster to the unit circle, the 
centroid (unweighted center of mass) is given in complex 
coordinates by the conjugate of the £ = 1 coefficient ijjl. 
Using a frequency term other than £ = 1 multiplies each 
angle ^ by a factor, effectively stretching or compress- 
ing the pattern along the circle. We see that choosing 
I = ipn-, where ipn is a rotational symmetry of the cluster, 
causes the different angles 9j to coincide, resulting in a 
non-vanising centroid for the transformed cluster. Thus, 
the Fourier coefficient with 1^1 represents the centroid 
of the stretched or compressed pattern. Although the 
position of the centroid is dependent on the cluster ori- 
entation, this distance from the origin to the centroid is 
invariant under rotations. Thus, we obtain a rotation- 
invariant descriptor by computing the magnitude of the 
coefficient. 

These properties of Fourier coeffcients have been ex- 
ploited in the context of bond order parameters [H [9]. 
For example, particular coefficient magnitudes, such as 
?/^6 (oi"^ analogously q4 and qg, see below) have been 
used directly as scalar order parameters [H [SQ. This 
method is often sufficient for the simple structures that 
we encounter in systems of spherical particles. However, 
more complex structures often require a full range of 
coefficients, and order parameters must be constructed 
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FIG. 4: Properties of harmonic descriptors, (a) The depicted 
cluster of four particles is somewhat 3-fold symmetric. When 
each angle is multiplied by a frequency factor = 3, the 
centroid position moves farther from the origin {\ip3\ > iV^il), 
as depicted by the red vectors. (6) When the cluster is ro- 
tated (i.e., a constant factor Oc is added to each angle 0), \ip3\ 
remains constant, although ips changes, (c) When the clus- 
ters on the left are scaled by a frequency factor ^ = 3, the 
transformed centroids ip3 become equivalent, even though the 
underlying shapes are different. 



by comparing sets of coefficients, rather than evaluating 
particular coefficients. To create a descriptor from the 
Fourier coefficients, we simply combine the desired ip£ or 
IV^^I into a long vector. For example, a general rotation- 
invariant shape descriptor that is applicable to patterns 
on the circle over a range of symmetries is given by: 



s^2=< |Vol,IV'il,-IV'^_J > 



(4) 



It is easy to imagine how different combinations of the 
Fourier descriptors can be used to create shape descrip- 
tors with different levels of robustness and sensitivity to 
particular symmetries. For many applications it is com- 
mon to include only coefficients with specific symmetries, 
or to use rotation-dependent coefficients. 

As outlined in the previous section, many 3d structures 
are well represented by patterns on the surface of the 
sphere. The analogy to the Id DFT on the 2d surface of 
the sphere is known as the discrete spherical harmonics 



transform (DSHT)[56], given by: 



-i 

j = 1,2, ...n^in (5) 



£=1 



The terms Yf^[Oj^(j)j] are spherical harmonics, defined 
by y^^((9,^) = NpP^" {cos 0)exp{im(l)). The term 
is a normalization factor ^/ {2i -\- — m)\ / {i -\- m) \ and 
is a Legendre polynomial. We see that the DSHT is 
similar to the DFT except for an additional term depend- 
ing on the polar angle 6. 

The Fourier coefficients for the DSHT are given by: 



^pts 



m 



1,...^ 



(6) 



Unlike the circular coefficients t/j^^ which are complex 
numbers, the spherical coefficients are 2^ + 1 dimen- 
sional complex vectors. Like the circular coefficients, the 
spherical Fourier coefficients are robust under noise, sen- 
sitive to rotational symmetries corresponding to ^, and 
can be used to construct invariants. Although the geo- 
metrical interpretation of these properties is more com- 
plex than for the Id case, the same principles apply. 

The rotation-invariant version of the spherical coeffi- 
cients is given by: 
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(7) 



Like the circular invariants the spherical invariants 
Iq^l are positive real numbers. 

In analogy with our example above, we can create a 
rotation-invariant shape descriptor for a pattern on the 
sphere by: 



S^^=< |qo|,|qi|,...|q,_J>. 



(8) 



Again, the optimal descriptor for a particular applica- 
tion depends on the desired properties such as robust- 
ness, sensitivity to specific symmetries, and rotation in- 
variance. 

Notice that although we use the notation and 
"q^" to highlight the connection with bond order param- 
eters, we have redefined the order parameters slightly, by 
changing the sign of the complex exponential (i.e., the 
conjugates in equations [2] and [6|. This sign change is 
inconsequential to the properties of the coefficients; our 
redefinition simply allows us to highlight the important 
relationship between tpi and q^ and the DFT, which is 
standard and extensively studied^]. An overview of the 
Fourier descriptors method is given in Fig. [4j As a gen- 
eral rule of thumb, is used when we only wish to 
describe the 2d ordering of a system (e.g., the spatial or- 
dering within a single plane in a confined fluid ^ |55j or 
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FIG. 5: Drawbacks of Fourier descriptors. Rot at ion- invariant 
Fourier descriptors for structures with multiple shells are in- 
sensitive to the relative orientation of inner and outer shells. 
Using invariant Fourier descriptors to describe the two struc- 
tures shown in the schematic shown erroneously produces 
identical results. 



the crystalline arrangement of cylindrical domain [4T, '58^) 
and "q^" is used when we wish to describe the 3d ordering 
(e.g., the spatial ordering in a 3d crystal [15^ or structure 
of a compact 3d cluster pT]). 



B. Zernike Descriptors 

As mentioned in the previous section, when patterns 
cannot be properly represented by the surface of a single 
circle or sphere, one solution is to break up the pattern 
into independent radial shells, and compute the Fourier 
descriptor for each shell independently. We then con- 
struct a shape descriptor by combining the Fourier de- 
scriptors for each shell into a long vector: 

SF, multishell ^ q^F q^F q^F ^ /q\ 

— ^ ^shelh^^shelh^ -"^shellr, ^ ' \^ ) 

Here, represents a Fourier descriptor, either S^^ or 
S^^, as defined in the previous section. 

While this scheme is sufficient for many problems, it 
has the drawback that small perturbations to the par- 
ticle positions can cause maxima in the pattern to shift 
between shells, causing errors in matching, particularly 
when Upts is small. Another drawback is that, since the 
shells are treated independently, the rot at ion- invariant 
Fourier descriptors are insensitive to relative orientations 
between the different shells; this is depicted in Fig. [5] for 
two structures that would erroneously have matching de- 
scriptors. 

We can solve these problems by representing our pat- 
tern in a coordinate system with radial dependence, 
which allows us to properly index patterns defined on 
the disk or ball rather than the unit circle or sphere. 
To do so, we use Zernike radial polynomials in our 
expansion [59l |6Q] . 

We can express the intensity of a pattern at a given 
point on the unit disk as: 

fi^j^^j) = ^^aniRni{rj)e:^-p [iiOj]. 

n i 

j = 1,2, . ..riMn (10) 



The terms Oj and rj represent the position of a point on 
the unit disk. The value of £ is restricted such that £ < n 
and (n — £) is an even number. The expansion coeffi- 
cients ani are known as "Zernike moments" and can be 
considered analogous to Fourier coefficients for the r, 6> 
coordinate system. The function Rniir) is a radial poly- 
nomial, where r is the radial distance from the center of 
the disk. Thus, the 2d Zernike expansion is very simi- 
lar to the Id Fourier expansion, but with an additional 
radial term. 

The Zernike moments are given by [59]: 

^ ^ ^ Upts 

ani = ^fj{rj,Oj)Rn£{rj)ex.p[-i£Oj]. (11) 

The terms Oj and Vj represent the position of an input 
point Xj in polar coordinates, normalized on the unit 
disk. Again, we require that £ < n and (n — £) is even. 
Each moment is a complex number. Since the radial 
polynomial is only dependent on r, the same invariance 
relations hold for the Zernike moments as for the Fourier 
coefficients. To define a rot at ion- invariant moment on 
the disk, we take the complex magnitude of the moment: 

\cinA = CLnldni- (12) 

The 2d Zernike invariants are positive real numbers. We 
can create a Zernike descriptor by combining many mo- 
ments in a vector. For example, we can create a 2d 
rotation-invariant Zernike descriptor by: 

S^^ =< |aoo|,|aiiUa2o|, 1^221,. ..|a^_.^_J > • (13) 

Like the Fourier descriptors, the frequency parameter 
£ has a straightforward relationship with the rotational 
symmetry of the pattern. Thus, we can sometimes choose 
important moments a priori. However, we typically take 
all moments with £ within a limiting frequency £max- 

We can express a pattern on the unit ball as the sum 
of 3d Zernike moments: 

/(^i, ^„ = E E E ZnMr,)Yr{9,,<i^j). 

n i m 

j = 1,2, ...riMn (14) 

The terms give the position of a point on the 

unit sphere. Again, we require £ < n and (n — £) is 
even. The moments are defined similarly to the Fourier 
coefficients on the surface of the sphere, but again with an 
additional radial component. The 3d Zernike moments 
are given by ^60^: 

m = -£,-£^ !,...£ (15) 

The variables [rj,^j,^j] represent the position of an in- 
put point Xj in spherical coordinates, normalized on the 
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unit sphere. Whereas the 2d Zernike moments a^^ are 
complex numbers and the 3d Zernike moments z^^ are 
complex vectors of length 2i -\- 1. Analogously to the 
spherical Fourier coefficients, we take the magnitude of 
the complex vector \zn£\ to define invariant moments on 
the unit ball: 



\Zn£\ 



2ZT1 



E K 



m 12 



(16) 



The 3d Zernike invariants are positive real numbers. Like 
the Fourier coefficients, they are sensitive to the rota- 
tional symmetries of the pattern and are robust under 
small perturbations. We can create a rotation invariant 
symmetry-independent 3d Zernike descriptor according 
to: 



= < |Zoo|, |Zii|, |Z2o|, |Z22|, ...|Z£ 



> . 



(17) 



When computing either multiple-shell Fourier descrip- 
tors or Zernike moments it is essential that the patterns 
being compared are normalized consistently. In the case 
of Zernike moments, all points in {x} must lie on the 
unit ball or disk. Typically, patterns are normalized by 
translating the centroid to the origin and rescaling the 
coordinates such that every point on the pattern has a 
radial distance less than 1. This scheme is sufficient for 
the majority of patterns that we encounter in particle 
systems. An overview of the Zernike scheme is depicted 
in Fig. [6| 
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C. Computational Considerations 

Fourier and Zernike coefficients can be computed from 
either point cloud data (i.e., raw particle positions) or 
voxel data (i.e., volumetric data or pixel data). As men- 
tioned previously, the two representations are essentially 
equivalent; point cloud data represent the limit of zero 
bin size, where each fi is equivalent ly 1. While this dis- 
tinction does not affect the properties or definition of 
the descriptors, it becomes important when considering 
the computational cost of a matching application. If the 
input data is point cloud data, we must compute the 
descriptors for each structure independently. However, 
for volumetric data, we can compute the contribution to 
each coefficient for each point on the grid beforehand, and 
then simply multiply by the intensity of each shape fi to 
compute the value of the coefficients. This can greatly 
reduce the computational cost when Upts is large or many 
coefficients are used. 

As an additional consideration, the time required for 
computing the transforms themselves can be greatly re- 
duced by computing the fast Fourier transform (FFT) 
rather than the DFT (or the equivalent for the appropri- 
ate coordinate system). Methods for computing the FFT 
and the discrete spherical harmonics transform, respec- 
tively, are given in references [57 and ^56^. An efficient 



FIG. 6: Zernike descriptors, (a) Decomposition of the pat- 
tern obtained for a 2d pentagonal cluster with Zernike mo- 
ments. Since the cluster has 5-fold symmetry, the leading 
coefficients are given by multiples of ^ = 6. (b) Decomposi- 
tion of the pattern obtained for a 3d fee cluster with Zernike 
moments. Since the cluster has 4 and 6-fold symmetry, the 
leading coefficients are given by multiples of ^ = 4 and i = 6. 
(c) Schematic of the resulting Zernike descriptors. The value 
of the each component in the descriptor is proportional to 
the contribution of the corresponding moment to the over- 
all sum. Notice that rotation-invariant Zernike descriptors 
contain real, positive components, while rotation-dependent 
descriptors contain complex components. 



method for computing Zernike coefficients is given in ref- 
erence [6Qj . 



III. QUANTIFYING SIMILARITY 

The shape descriptors derived in the previous section 
can be considered compact mathematical representations 
of the underlying particle structures. The physical sim- 
ilarity between different particle structures can then be 
quantified by the mathematical similarity between shape 
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descriptors. Shape descriptor similarity is quantified by 
a similarity metric "M," that gives a scalar value that is 
proportional to the similarity between descriptor pairs. 
For convenience, we define M such that, by construc- 
tion, it lies on the interval M G [0,1], or sometimes 
M G [—1,1]. This strengthens the analogy between M 
and what we normally consider to be an order parameter, 
since order parameters typically give a value of 1 for per- 
fectly ordered structures and for perfectly disordered 
structures. Since our harmonic shape descriptors are de- 
fined as vectors, we can define similarity metrics based 
on standard vector operations. For example, one sim- 
ple similarity metric is given by the Euclidean distance 
between shape descriptor vectors: 

M,,,,(Si, S2) = 1 - |Si - S2|/(|Si| + IS2I). (18) 

Variations of Mdist are common throughout the shape 
matching literature. Another simple similarity metric is 
proportional to the dot product between shape descrip- 
tors: 

M,,,(Si,S2) = Si-S2/(|Si||S2|). (19) 

Schemes similar to Mdot have been used in applications 
involving standard bond order parameters, such as mea- 
suring correlation lengths [9l |T2] and identifying crystal 
grains [15 . 

Similar information is obtained from Mdist and Mdou 
the only difference is that whereas Mdist is more sensitive 
to the absolute difference between vector components, 
Mdot is more sensitive to the overall direction of the vec- 
tor and the sign of the components. Thus, Mdot is often 
superior when matching non-ideal structures from a par- 
ticle system to mathematically perfect reference struc- 
tures, since thermal noise will tend to damp the frequency 
domain signal, but the descriptor will retain the same 
basic character for a given class of structures and hence 
the same direction. The Mdot metric may also be favor- 
able when comparing rotation-dependent harmonic de- 
scriptors, which may contain either negative or positive 
components, whereas Mdist iiiay be favorable for invari- 
ant descriptors, where all values are inherently positive. 
Since compared shapes are usually at least grossly simi- 
lar, matching values are rarely for either metric. As is 
discussed in greater detail in reference [35 , it is often nec- 
essary to determine a lower bound on M by comparing 
to structures that are known to match poorly to obtain 
a baseline value. 



IV. EXAMPLE APPLICATIONS 

Shape similarity information obtained from evaluating 
the match, M, between shape descriptors, can be applied 
to create various types of structural metrics for complex 
particle systems. In this section, we provide several ex- 
ample applications that can be addressed by the use of 
shape matching with harmonic descriptors. We provide 




FIG. 7: Tracking structural transitions in the protein Ubiq- 
uitin as it is pulled from both ends 62| |63]. The protein is 
pulled linearly as a function of time. The black curve shows 
a matching order parameter using the initial folded configu- 
ration i as a reference structure. The red curve shows the 
order parameter using the final unfolded configuration / as 
the reference structure. The blue dashed lines highlight the 
times at which significant structural changes occur. 



a more extensive range of example applications based on 
shape matching methods in general in reference [35 . Sev- 
eral additional examples of shape matching applications 
based on the descriptor are given in reference [34 . 



A. Order Parameters and Correlation Functions 

Perhaps the most standard application of order pa- 
rameters is to track structural transitions, either as a 
function of time or a changing reaction coordinate. As 
an example, consider the protein "Ubiquitin[61 " shown 
in Fig. [7| which unfolds as it is pulled from both ends. 
This is a standard example problem from the NAMD and 
VMD tutorials [621 ES]: both of which are available on- 
line. Since the structure is 3-dimensional and has radial 
dependence, we can index it using a Zernike descriptor 
on the unit ball, S^^. To match the shape independently 
of the orientation of the sheet, we take rotation invariant 
moments with I in the range 4 < i < 12. We use both 
the initial folded state i and the final unfolded state / as 
reference states. Fig. [7| shows the unfolding transition as 
a function of time tin a NAMD molecular dynamics sim- 
ulation. The protein unfolds in three steps (blue dashed 
lines), in agreement with visual inspection. The noise in 
the data is indicative of the thermal fiuctuations in shape 
observed at this temperature. 

As a slightly more complex example of characteriz- 
ing transitions, consider the micro-phase separated struc- 
tures formed by ditethered nanospheres [58] shown in 
Fig. |8] The structure of the equilibrium system goes 
through two transitions as a function of the effective in- 
verse temperature, first from a disordered structure to 
a tetragonal cylinder/tetragonal- mesh (TC/TM) phase 



10 





1 1 1 
M(S, Sdis) 


1 1 1 1 1 1 1 1 
M(S, Stc/tm) I M(S, Stc/tc) 


1 — *i 


1 1 1 1 1 ^ \ \ — 



0.25 0.5 0.75 1 1.25 1.5 
l/keT (strength of attraction) 

FIG. 8: Transitions in ditethered nanospheres. (a) Depiction 
of the three structures formed by a ditethered nanosphere 
system [58] (left to right: disordered, TC/TM, TC/TC). No- 
tice that the cyUndrical structures span the simulation box in 
the z-dimension, into the page. (6) Matching order parame- 
ter for the three reference structures as a function of inverse 
temperature. 



and then to a similar tetragonal cylinder (TC/TC) phase. 
The abbreviations indicate the structure of the tethers 
(blue, red) and nanoparticles (white), respectively. We 
can quantify this behavior by matching the global pat- 
terns obtained at different temperatures with ideal struc- 
tures taken from within the three structural regimes: the 
disordered regime, the TC/TM regime, and the TC/TC 
regime. The global pattern for each structure is char- 
acterized by the probability distribution of local density 
maps for each particle type, as depicted in Fig. |2]i. 
To capture ordering on a range of length scales, den- 
sity maps are computed for four radial shells with r = 
[3cr, 4.5cr, . . . 9cr], where a is the characteristic lengthscale, 
given by the diameter of the tether beads in the simula- 
tion. For each shell, we compute the rot at ion- invariant 
Fourier descriptors S^^, where we take a range of frequen- 
cies 4 < ^ < 12. A pseudo-order parameter for each refer- 
ence structure is then given by M^if f{S^^ ^S^^f). Fig. [7| 
shows the order parameters for the three reference struc- 
tures as a function of inverse temperature. We observe 
that the structural transition between the three phases 
is smooth and continuous, as verified by visual inspec- 
tion. In reference [35 , we show that, for this particular 
problem, we can obtain a nearly identical result using a 
simpler shape descriptor akin to the radial distribution 
function g{r). However, for more complex phase sepa- 
rated structures, harmonic descriptors typically give a 
better representation of the underlying shapes than such 
coarse measures as g{r). 

In addition to characterizing how structures change 
as a function of time or a reaction coordinate, another 
common application is to characterize how structures 



change in space by computing correlation functions. In 
this case, we choose structures from different points in 
space, rather than ideal structures, as references. Sev- 
eral examples of spatial correlation functions based on 
bond order parameters have already been defined in 
the context of measuring lengthscales for crystal-like 
ordering [8^, [Q] |14| |64]. This involves measuring quanti- 
ties such as (M(Si, Sj)) (nj), that give the average sim- 
ilarity value as a function of the radial separation rij. 
Typically, a particular rotation-dependent Fourier coef- 
ficient, for example ipQ or q^, is chosen for S and M^ot 
is chosen for M. An alternative, and well known, spatial 
correlation function based on bond order parameters is 
given by the • crystal grain detection scheme of refer- 
ence [15 , which has been applied to studying nucleation 
and growth and characterizing crystalline defects in sev- 
eral simulation and experimental studies involving close- 
packed and bcc crystals[16, ilTj l20J ES]. Here, crystal 
grains are identified within a bulk liquid by first notic- 
ing that, for many crystals, local clusters within crys- 
tal grains match with their neighbors in terms of both 
shape and orientation, whereas clusters in the liquid do 
not. Thus, pairs of particles j in grains typically satisfy 
M(Si,Sj) > Mcut^ where S is a lotdition- dependent har- 
monic descriptor and Mcut is a sufficiently high value so 
as to exclude poor matches. Even in the liquid, some 
pairs of particles inevitably satisfy M(Si,Sj) > Mcut 



due to random fluctuations. Thus, 
crystal-like ordering is given by: 



a local indicator of 



1 



rinbr 

E 



e[M{Si,Sj) - Me. 



(20) 



Here, O is the Heaviside function. Typically, we enforce 
fi > fcut , where the value of fcut is chosen to distinguish 
liquid and crystal-like configurations [1 5, T^. Although 
the original scheme is based on the Fourier coefficient 
as a shape descriptor (for identifying fee, hep, and bcc 
crystals) and the M^ot similarity metric, other instances 
of S and M can be used depending on the structure under 
investigation. For example, in references [43j and [19], we 
used different bond order parameters as shape descrip- 
tors to identify particles in the diamond lattice and in 
dodecagonal quasicrystals, respectively. Arbitrarily com- 
plex crystal structures can be treated using this method 
by harmonic descriptors with a full spectrum of Fourier 
coefficients or Zernike moments. We apply this scheme 
in the context of two examples highlighting the special 
symmetry properties of harmonic descriptors in section[V| 
below. 



B. Database Search and Structure Identification 

In addition to computing order parameters and cor- 
relation functions, another common application of bond 
order parameters is to identify local structures such as 
icosahedral clusters [lOl [TTJ [191 [41] within a bulk system. 
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FIG. 9: Structure identification using a tunneling electron 
microscopy (TEM) image from reference [5]. The structure 
depicted in the image is identified by finding the best match 
from a database of candidate structures. For this particular 
example, the reference structures are created by mathematical 
construction. The image is correctly identified as depicting 
the CuAu structure viewed along the (100) direction. 



This typically involves choosing a cutoff value for a par- 
ticular Fourier coefficient (for example, l^el) above which 
a cluster is identified as the structure of interest. This 
is, in its essence, a rudimentary shape matching scheme, 
where the Fourier coefficient provides a coarse description 
of the cluster shape, and the cutoflF acts as a similarity 
metric. This type of structure identification scheme can 
be applied within a much broader context by using har- 
monic descriptors to perform a database search for an 
unknown structure. The unknown structure is identified 
as the structure from the database of known structures 
that gives the best match. Database searches based on 
harmonic descriptors have already been applied to pro- 
teins and macromolecules[36^i |37]. In an earlier publi- 
cation, we applied a database search to identify local 
structures within a phase separated system of tethered 
nanoparticles[ll . In the future, they may be applied to 
more abstract problems, such as data mining for web- 
accessible particle structures. 

Database searches based on harmonic descriptors can 
be applied to a wide range of complex local and global 
structures. As an example, consider the tunneling elec- 
tron microscopy (TEM) image depicted in Fig. [9j which 
shows nanoparticles arranged in a binary crystalline su- 
perlattice from reference [5 . Although this structure was 
identified as the AuCu crystal viewed along the (100) di- 
rection by visual inspection [5 , assume for the purpose of 
this example that the structure is unknown. The struc- 
ture of the lattice can be identified by finding a best 
match from a database containing the images of known 
reference structures. For our proof-of-concept example. 



we use a minimal reference database consisting of four 
different ideal binary crystal structures; however for more 
realistic problems, the database may be much more ex- 
pansive. The reference structures are created by mathe- 
matical construction and rendered by placing spheres at 
the lattice positions. In practice, matching can be per- 
formed using other non-ideal images or other experimen- 
tal images. The images are indexed for comparison using 
the 2d Zernike descriptor S^^. As mentioned previously, 
harmonic descriptors can be used to describe images [59]. 
where {x} and {/} represent the positions and intensities 
of individual pixels. In practice, the pixel intensities are 
inverted (/^ = — /i,o), since, for the current set of images, 
the particles are darker than the background. To ensure 
that the matching algorithm is not affected by the differ- 
ent particle shades, we apply a binary thresholding crite- 
rion li = Q{Ii — I cut)' To extract a global pattern from 
the image, we use the probability distributions method 
depicted in Fig. [l] (Notice that although the structures 
are crystalline, the superposition method is not appli- 
cable, since the particle centers are not known). For 
each local structure, we compute S^^ descriptors with 
rotation-invariant moments and frequencies in the range 
< ^ < 8. Our overall results are not impacted by 
the inclusion of higher frequencies. For each image, lo- 
cal descriptors are computed for 100 different randomly 
chosen pixels. For each pixel, the range of neighboring 
pixels used to construct the local descriptor corresponds 
to roughly three particle diameters. The local descrip- 
tors are then combined into an overall probability his- 
togram descriptor for matching. As shown in Fig. |9j the 
unknown structure most closely resembles the CuAu lat- 
tice along the (100) direction, in agreement with visual 
inspect ion [5]. Additional orientations of the crystalline 
lattices could also be considered to better identify struc- 
tures with complex ordering. 

This same basic identification scheme can be used for 
all types of structures, either simulated or experimental, 
local or global. A subtlety arises when disordered local 
structures are possible, since it is typically infeasible to 
construct a reference database for the vast space of "dis- 
ordered" structures. As outlined in references [TTJ [35] 
this problem can be avoided by setting a minimum value 
for the best match, below which structures are considered 
disordered. 

In addition to database searches, experimental images 
can be used for all of the other applications described here 
or in references [34 and [35 , including quantifying struc- 
tural perfection, computing order parameters and corre- 
lation functions, grouping and classifying structures, etc. 
The only caveat is that, like simulation data, the images 
must be properly normalized such that the particle sizes 
and lattice spacings (if applicable) are the same for all 
compared structures. In addition to images of particles, 
information can be extracted from other types of images, 
such as diffraction patterns. For many applications, im- 
age processing can be applied to experimental images 
to simplify the data j66l - [68] . for example, by identifying 



12 



particle centroids. We explore the combination of image 
processing techniques and shape matching algorithms in 
a separate publication [69]. 



V. SPECIAL PROPERTIES OF HARMONIC 
DESCRIPTORS 



In the previous section, we applied harmonic descrip- 
tors within the context of general particle shape matching 
applications, similar to those outlined in reference [35^. 
Thus, although the harmonic descriptors have useful 
properties, such as rotation invariance, we could just as 
easily base our examples on other shape descriptors with 
similar properties. In contrast, in this section we ex- 
plore applications for which harmonic descriptors and 
Fourier coefficients are specifically well-suited, due to 
their unique symmetry properties. 



A. Matching to Within an n-Fold Rotation 

One such application is the problem of matching struc- 
tures that are a unique rotation of one another. As an ex- 
ample, consider the problem of detecting the local crystal 
grains of different structures in the 2d system depicted in 
Fig. [TO^, where particles interact via the Lennard- Jones 
Gauss (LJG) potential [70]. For a potential-minimum- 
position parameter ro = 1.57, the system forms two 
crystal structures: the honeycomb (he) structure at high 
values of the well-depth parameter e and the hexagonal 



(hex) structure at low e[70 . As outlined in section IV A 



local crystal-like pairs are typically identified by local 
clusters that match in terms of both shape and orienta- 
tion. This method is easily applicable to the hex lattice; 
however, in the case of the he lattice, the triangular first 
neighbor shells of neighboring particles are mirrored and 
rotated by 60° in the plane (see Fig. [ToJd). Thus, they 
match in terms of shape, but not in terms of orienta- 
tion. The symmetry properties of the Fourier coefficients 
pose a unique solution to this problem. In the space 
of the ^ = 3 Fourier coefficient ?/^3, the triangular neig- 
bor shells in the he lattice are precisely antiparallel (i.e., 
Mdot[^2.(i), V^(j)3] = -1, where i and j are neighbors). 
Thus, a matching criterion can be constructed based on 
V^s) to determine whether two neighbor shells 



are in the ideal he configuration. In Fig. 10 i, particles 
in the he structure are colored blue, particles in the hex 
structure are colored red, and other particles, that do not 
belong to a particular crystal grain, are colored black. 
The 3d analogy of this method, using in the place of 
?/^3, was used to measure the number of particles in local 
diamond lattice grains in a system of patchy particles in 
reference [43 (see Fig. flOp). 




FIG. 10: Honeycomb (he) to hexagonal (hex) transition in the 
2d Lennard- Jones Gauss (LJG) system [70 . (a) Annealing re- 
sults obtained for a LJG well- minimum position ro = 1.57 for 
three different values of the well depth parameter s. From top 
to bottom: e = 2.7, 2.2, 1.4. Particles in the he structure are 
colored blue, while particles in the hex structure are colored 
red. All other particles are colored black. (6) Depiction of the 
local neighbor configurations in the hex and he structures, (c) 
Diamond lattice formed by patchy particles [431 171]. Analo- 
gously to the he lattice, the diamond lattice can be character- 
ized by the negated match of the third-order harmonic [43] . 



B. Matching Based on Rotational Symmetries 

A similar application for which Fourier coeffcients are 
uniquely suited is the problem of matching structures 
based on their rotational symmetries rather than their 
shapes. As an example, consider the problem of match- 
ing local neighbor shells in the decagonal (i.e., 10-fold 
symmetric) quasicrystal formed in the 2d LJG system [TOj 
(Fig.jn^). Over the range indicated, the neighbor shells 
exhibit strong 10- fold rotational symmetry with a com- 
mon direction, but the neighbor shells have different 
shapes. Thus, our criterion for detecting local crystal 
grains outlined in section IV A| fails for most shape de- 
scriptors, since the underlying shapes do not match. As 
a solution, we can describe each local cluster with the 
^ = 10 Fourier coefficient i/jiq. Since the clusters are 10- 
fold symmetric, and oriented in the same direction, the 
complex number iJjiq is identical regardless of whether 
the clusters are missing particles. Local quasicrystalline 
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grains can then be detected as outlined in section |IV A[ 
using Mdot {i^io^ '010 ) to identify local crystal-like pairs. 

In reference [19] we use an analogous scheme to detect 
ordered grains in a 3d dodecagonal (12-fold symmetric) 
quasicrystal (Fig. p!T]3). In this case, the structure has 
hundreds of different neighbor shell directions [42], which 
exhibit strong 12-fold symmetry. This is depicted by the 
superposition of local patterns over the range r < 2. 3 la 
(i.e., ~ 2 neighbor shells) in Fig. 11 z. Considering this 



longer range ensures that each local region contains a 
sizable fraction of 12-fold directions. Since we are only 
interested in the rotational symmetry and directional- 
ity of these local patterns, we remove the r-dependence 
from the patterns prior to matching (Fig. [Tl]^). To cap- 
ture the 12-fold symmetry, we use a matching criterion 
based on Mdot{^i2^ ^12)- Particles with a minimal frac- 
tion of solid-like matches fcut are considered to be lo- 
cally quasicrystalline. The cutoffs are determined by tak- 
ing the crossover points in the probability distributions 
^(Mrfot(qi2,qi2)) and P{f solid) (Fig. 
reference [T9^, we take Mcut = 0.45, fcut 



Hi). Following 
0.5. As de- 



picted in Fig.pT^, this criterion is sufficient to determine 
a small quasicrystal nucleus in the bulk liquid. 

Notice that although both of our examples here are 
based on quasicrystalline structures, the method of 
matching dissimilar structures based on their rotational 
symmetries is applicable to a wide range of structures. 
As a trivial example, we return to our previous prob- 
lem of detecting crystal structures in the 2d LJG system. 
Suppose now that we require an order parameter that 
detects crystalline grains in general, either hex or he. 
In this case, we can use a scheme based on the i = 6 
Fourier descriptor 7/^65 since in the space of ^ = 6 har- 
monics, the triangle neighbor shells of the he structure 
and the hexagon neighbor shells of the hex structure are 
equivalent. 



C. Orientation About a Symmetry Axis 
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FIG. 11: Matching based on rotational symmetries, (a) 
Decagonal (10-fold) quasicrystal formed in the 2d LJG 
system [70]. The call-out depicts local structures within the 
system, which are 10-fold symmetric but dissimilar in shape, 
(a) Dodecagonal (12-fold) quasicrystal formed in the 3d 
Dzugutov system by MD simulation [72 , with a simulated 
diffraction image (lower right), (b) superposition of all neigh- 
bor clusters with a cutoff range 2.31cr. The images on the 
right disregard r-dependence. (c) Probability distributions of 
bond correlations and solid-like neighbors used to determine 
the criterion for local quasicrystal grains, (d) Dzugutov liquid 
with a single small quasicrystal nucleus of about 50 atoms. 



As a final example of an application that exploits the 
symmetry properties of Fourier coefficients, consider the 
problem of aligning a crystal structure about a particular 
symmetry axis. This is useful, for example, for computer 
algorithms that depend on the orientational direction of 
a structure, such as automatically computing the diffrac- 
tion pattern or rendering images of particle data. Sup- 
pose, for example that the desired symmetry axis is the 
i = io plane, and the desired alignment direction is the 
z-axis. The crystal can be iteratively rotated to maxi- 
mize the Fourier coefficient \^piQ\ where all particles are 
projected onto the xy plane. If the crystal consists of a 
single grain, a relatively small test cluster can be used 
to perform the optimization, greatly enhancing compu- 
tational efficiency. 

As an example, consider the problem of aligning the 
face-centered cubic (fee) crystal shown in Fig. 12 1 along 
an 8-fold planar symmetry axis. Finding the optimal 



rotation that maximizes lipg] in the plane is solved by 
applying a simple simulated annealing Metropolis Monte 
Carlo (MC) scheme. Our MC scheme involves attempt- 
ing trial rotations of a small test cluster from the cen- 
ter of the box, which are accepted according to a Boltz- 
mann probability distribution exp[— /3|?/^8|], with the fic- 
titious energy function — IV^sl- The inverse temperature 
P is increased from 0.5 to 1000 over 100 steps. In prac- 
tice, the scheme only converges to a local minimum in 
the energy — IV^sl, so many different initial orientations 
are attempted to find a global minimum. As depicted 
in Fig. [T2j optimizing lips] aligns the 8- fold symmetry 
axis of the structure along the z-axis, but the structure 
is misaligned with the x and y-axes. To perform this 
alignment, we zero the complex component ^(^/^g)^ such 
that 5R('08) = IV^sl- As depicted in Fig. 12 3, this scheme 
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FIG. 12: Orientation about a symmetry axis, (a) Depiction of an algorithm for orienting a face- centered cubic (fee) crystal 
about an 8-fold planar symmetry axis, in this case, for the purpose of automatically computing the diffraction image. First, 
the planar symmetry, measured by \ip8\, is maximized such that the 8-fold symmetry axis is oriented in the z-direciton (out of 
the plane). Then the structure is rotated in the xy plane to zero the complex component ^(ipg), such that the symmetry axis 
is oriented with the xy axis. Finally, the diffraction image is computed, (b) The same process, repeated for a structure with 
thermal disorder. 



is robust, even under a large amount of thermal noise. 
More complex optimization algorithms may be applied to 
improve computational efficiency and accuracy over our 
simple MC scheme. This type of orientation algorithm 
is potentially useful for matching diffraction data, since 
a large number of simulated diffraction images for 3d 
structures can easily be computed automatically about 
various symmetry axes. 

VI. SUMMARY AND FUTURE OUTLOOK 

In summary, we have demonstrated how bond order 
parameters, already defined for particle structures on the 
unit circle and sphere, can be extended to index struc- 
tures on the unit disk or ball. We have demonstrated how 
these bond order parameters can be used to create har- 
monic shape descriptors, which can in turn be applied to 
create unique, highly specialized order parameters and 
automatically identify unknown particle structures. In 
addition to the minimal proofs of concept reviewed here, 
more complex matching applications are explored in ref- 
erence [35] . 

In the future, the ability to describe structures nu- 
merically lends itself to many novel applications. In the 
short term, matching applications can be used to auto- 
mate structural analysis for large datasets. Shape de- 
scriptors can be used within the context of many of the 
enhanced computational algorithms used in self-assembly 



and computational biology, such as path sampling [7?. [74] 
or metadynamics[75 in the context of pseudo order pa- 
rameters or collective variables to guide the sampling. 
Shape descriptors can also serve as the basis for new 
optimization algorithms such as genetic algorithms [76j, 
which often rely on energy rather than structure as a nu- 
merical measure of fitness. In addition to abstract com- 
putational applications, shape matching algorithms can 
be applied to experimental images to obtain quantita- 
tive insight into experimental data. By combining shape 
matching algorithms with new image processing schemes, 
much of the same information that we have obtained for 
simulation data can be obtained for experimental systems 
as well. 
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